function val = tension_spline(a, x, y)
  % @param
  %   x = vectorul absciselor punctelor din suportul de interpolare
  %   y = vectorul ordonatelor punctelor din suportul de interpolare
  %   a = vectorul punctelor in care vreau sa calculez valorile functiei
  % @return
  %   val = vectorul ce contine valorile functiei in punctele din a
  [m n] = size(x);
  val = zeros(1, length(a));
  h = diff(x);

  %construiesc matricea
  A = diag(h(1:n-1), 1);
  A = A + diag(h(1:n-1), -1);
  A = A + 2 * diag([h(1), h(1:n-2) + h(2:n-1), h(1)]);

  %construiesc vectorul termenilor liberi
  b = 3 * diff(y(1:n));
  b = b / h(2);
  b(1) = b(1) * h(2) / h(1);
  b = [b,3 * (y(n)-y(n-1)) / (x(n)-x(n-1))];
  aux = 3 * diff(y(1:n));
  aux = aux / h(1);
  aux(n-1) = aux(n-1) * h(1) / h(n-1);
  aux = [3 * (y(2)-y(1)) / (x(2)-x(1)), aux];
  b = b - aux;
  b = b';

  %rezolv sistemul obtinut => vectorul c
  c = A\b;

  %construiesc vectorul b
  b = diff(y(1:n));
  b = b ./ h(1:n-1);
  aux = zeros(1, n-1);
  aux(1:n-1) = (h(1:n-1))' / 3 .* (2 * c(1:n-1) + c(2:n));
  b = b - aux;

  %construiesc vectorul d
  d = diff(c(1:n));
  d = d' ./ (3 * h(1:n-1));

  %valoarea in a
  for i=1:length(a)
    aux = abs(x(2:n) - a(i));
    %ind = indicele minimului
    [xmin ind] = min(aux);
    if a(i) > x(ind+1) && a(i) < x(ind+2)
      ind = ind + 1;
    endif
    val(i) = y(ind) + b(ind)*(a(i)-x(ind))+c(ind)*(a(i)-x(ind))^2+d(ind)*(a(i)-x(ind))^3;
  endfor
endfunction
